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Abstract 

In this paper a new technique aimed to obtain accurate estimates of the 
error in energy norm using a moving least squares (MLS) recovery-based pro- 
cedure is presented. We explore the capabilities of a recovery technique based 
on an enhanced MLS fitting, which directly provides continuous interpolated 
fields, to obtain estimates of the error in energy norm as an alternative to 
the super convergent patch recovery (SPR). Boundary equilibrium is enforced 
using a nearest point approach that modifies the MLS functional. Lagrange 
multipliers are used to impose a nearly exact satisfaction of the internal equi- 
librium equation. The numerical results show the high accuracy of the pro- 
posed error estimator. 

KEYWORDS: error estimation; equilibrated stresses; stress recovery; extended finite ele- 
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1 Introduction 

During the last few decades numerical techniques, such as the finite element method 
(FEM), have been used to approximate the solution of real problems. In order to 
assess the quality of these approximations it is necessary to evaluate the error ob- 
tained in the simulation. Current methods used to estimate the discretization error 
of finite element (FE) solutions are usually classified into different families: residual 
based, recovery based, dual analysis techniques,... [1, 2, 3]. The use of recovery 
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based estimators is widespread due to their robustness and simple implementation 
into existing FE codes. 

Today, novel techniques such as the extended finite element method (XFEM) are 
being used to introduce a priori known information about the problem solution into 
the FE formulation. The XFEM [4] enriches the classical FEM basis functions using 
a partition of unity approach in order to capture the local features of the solution 
in a cracked domain, i.e. the discontinuity of the displacement field along the crack 
faces and the singularity of the stresses in the vicinity of the crack tip. 

Although the XFEM provides highly accurate solutions, and significantly im- 
proves the modelling of certain types of problems, there is an urge to develop error 
control techniques for these kind of methods, mostly because of their increasing 
importance and the fact that they use rather coarse discretizations. For example, 
recovery based error estimators for partition of unity methods have been developed 
in [5, 6, 7, 8], using the residual approach in [9, 10] and the constitutive relation 
error (CRE) in [11]. In [12] the CRE is used for goal oriented error estimation in 
XFEM. 

The enforcement of the internal and the boundary equilibrium equations for 
stress recovery has been previously considered in the literature as a mean to im- 
prove the quality of the recovered field. For patch based formulations, [13, 14] 
introduced the squares of the residuals of the equilibrium equations to the least 
squares functional solved in the recovery process through a penalty parameter. In 
[15] a point-wise enforcement of the internal equilibrium in the polynomial basis, 
used to represent the recovered stress at the support of each node, was presented. 
Then, boundary equilibrium conditions were applied on a set of sampling points 
in the part of the patch boundary that coincides with the domain boundary. The 
SPR-C technique proposed in [16] imposed equilibrium constraints to the polyno- 
mial basis via Lagrange multipliers. The internal equilibrium was exactly satisfied 
at each patch, and a Taylor expansion of the applied stresses was enforced along 
the Neumann contour. Then, a conjoint polynomial procedure [14] was used to ob- 
tain a continuous stress field. Later, in [7] this technique was extended to XFEM 
approximations. 

Procedures to smooth or to recover the stress field based on MLS have also been 
used. In [17] a continuous stress field was obtained through local interpolation of the 
nodal displacement values using MLS. In [5] this same formulation was extended to 
XFEM problems, considering an enriched MLS basis and a diffraction criterion, and 
an error estimate for enriched approximations was proposed. In [18] a procedure 
to smooth the stresses for the meshless element free Galerkin (EFG) method using 
MLS shape functions was described. In [19] a so called Statically Admissible Stress 
Recovery Technique (SAR) that used MLS to fit the stress at sampling points was 
presented. The SAR technique comprised basis functions which consider the in- 
ternal equilibrium equations and the local tractions conditions along the Neumann 
boundary. In [20] an extension of the SAR technique for XFEM was presented. 
Following the definition of pseudo-divergence-free field used in [21], we can say that 
in references [19] and [20] the authors considered a pseudo- satisfaction of the in- 
ternal equilibrium equation. They indicated that accurate stresses were obtained 
with SAR, but they did not go further to evaluate any error estimate. In [22] dual 
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techniques were used in meshless methods to obtain an equilibrated dual problem 
using MLS shape functions that approximate Airy stress functions. 

The objective of this paper is to present an enhanced version of the MLS recovery 
technique to evaluate accurate estimates of the discretization error for FEM and 
XFEM problems that is based in the ideas presented in [7, 8, 16, 23]. The rationale 
behind the proposed technique is to try to enforce the recovered stress field to 
satisfy continuity (this property is directly provided following the MLS approach) 
and the equilibrium equations that are satisfied by the exact solution such that 
recovered stresses get closer to the exact stress field. An appropriate application 
of the equilibrium constraints is required to avoid discontinuities in the recovered 
field. For that reason, a novel approach to introduce the internal and boundary 
equilibrium constraints is proposed. The procedure has been implemented in a FE 
code where mesh refinement is based on element splitting and the use of constrain 
equations (Multi Point Constraints, MPC) to force C° continuity at hanging nodes. 
SPR requires special treatment of these nodes because it is based on the mesh 
topology. The use of the proposed technique is more flexible as it is not constrained 
by the topology of the finite element mesh. This feature is very powerful and it 
allows the direct use of the technique with isogeometric analysis with /i-adaptive 
refinement based on T-splines and in cases where the FE mesh is missing, like with 
meshless methods or elements with an arbitrary number of sides. 

Reference [23] showed that upper bounds of the error in energy norm can be 
obtained with recovery-based error estimators if the recovered stress field is statically 
admissible. The recovered stresses resulting from the use of the technique proposed 
in this paper are continuous, satisfy the contour equilibrium equation and provide 
a nearly exact satisfaction of the internal equilibrium equation (more accurate than 
the pseudo-satisfaction of the equilibrium equation used in previous works [19, 20]). 
Although the upper bound property is not guaranteed, the numerical results show 
that the proposed technique yields sharp error estimates which nearly bound the 
exact error. 

The paper is organised as follows: in Section 2 we present the reference prob- 
lems and their approximate solutions using FEM and XFEM. Section 3 deals with 
the main aspects of error estimation in FE approximations and the moving least 
squares formulation considering equilibrium conditions. Finally, numerical results 
are presented in Section 4, and conclusions are drawn in Section 5. 

2 Problem Statement and Solution 
2.1 Problem statement 

Let us consider the 2D linear elasticity problem. The unknown displacement field 
u, taking values in Q C M 2 , is the solution of the boundary value problem given by 



V • a (u) = b 
a (u) • n = t 
u = 



onT 



in Q 



on r 



D 



N 
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(2) 
(3) 



3 



where T N and T D denote the Neumann and Dirichlet boundaries with dfl = T N UT D 
and r^y H = 0. The Dirichlet boundary condition in (3) is taken homogeneous 
for the sake of simplicity. 

The weak form of the problem reads: Find u £ V such that 

VveV a(u,v) = Z(v), (4) 

where V is the standard test space for the elasticity problem such that V — {v | v £ 
^ 1 (^),v| rD (x) = 0}, and 

a(u,v):= / a T (u)e(v)dn= [ o-(u) T D-V(v)dft (5) 

Z(v) := / b T vdn+ [ t T vrfr, (6) 

where a and e denote the stresses and strains, and D is the elasticity matrix of the 
constitutive relation a = Ds. 

2.1.1 Singular problem: 

Figure 1 shows a portion of an elastic body with a reentrant corner (or V- notch), 
subjected to tractions on remote boundaries. For this kind of problems, the stress 
field exhibits a singular behaviour at the notch vertex. 




Figure 1: Sharp reentrant corner in an infinite half-space. 

The analytical solution of the stress distribution in the vicinity of the singular 
point is a linear combination of singular and non-singular terms. It is often claimed 
that the term with a highest order of singularity dominates over the rest of terms 
in a sufficiently close zone surrounding the singular point. The analytical solution 
to this singular elastic problem in the vicinity of the singular point can be found in 
[24, 25]. If a — 2tt the problem corresponds to the classic crack problem of linear 
elastic fracture mechanics (LEFM) that will be used in the numerical examples. The 
following expressions show the first term of the asymptotic expansion of the solution 
for mixed mode loading conditions in a 2D cracked domain [25]: 
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where and K\\ are the generalise stress intensity factors (GSIF) for modes I 
and II. The GSIF are multiplicative constants that depend on the loading of the 
problem and linearly determine the intensity of the displacement and stress fields in 
the vicinity of the singular point. 



2.2 Solution with FEM/XFEM. 

Let be a finite element approximation to u such that u^(x) = X^e/^( x ) u *> 
where N { represent the shape functions associated with node i and I is the set of all 
the nodes in the mesh. The solution lies in a functional space V h C V associated 
with a mesh of isoparametric finite elements of characteristic size h, and it is such 
that 

VvG^ a(iAv) = /(v) (9) 

Considering an XFEM formulation for the case of the singular problems of LEFM 
above-mentioned, the FE approximation is enriched with Heaviside functions to 
describe the discontinuity of the displacement field and with crack tip functions to 
represent the asymptotic behaviour of the stress field near the crack tip. This avoids 
the need for a conforming mesh to describe the geometry of the crack [26] and the 
use of adaptive techniques in order to capture the special features of the solution. To 
ensure the continuity of the solution, the partition of unity property of the classical 
linear shape functions is used. Therefore, the XFEM displacements interpolation in 
a 2D model is given by: 

u fc (x) = Y, N i(*)*i + E NjWHWbj + Y N ™(x) ( E F <M<) (10) 

iei jeJ meM \£=1 J 

where are the conventional nodal degrees of freedom, bj are the coefficients as- 
sociated with the discontinuous enrichment functions, and c m those associated with 
the functions spanning the asymptotic field. In the above equation, I is the set of 
all the nodes in the mesh, M is the subset of nodes enriched with crack tip func- 
tions, and J is the subset of nodes enriched with the discontinuous enrichment (see 
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Figure 2). In (10) the Heaviside function H ) with unitary modulus and a change of 
sign on the crack face, describes the displacement discontinuity if the finite element 
is intersected by the crack. are the set of branch functions used to represent the 
asymptotic expansion of the displacement field around the crack tip seen in (7). The 
Fa functions used in this paper for the 2D case are [26] : 




Figure 2: Classification of nodes in XFEM. Fixed enrichment area of radius r e 

More details regarding the XFEM implementation are given in the numerical 
examples section. 



3 Nearly Equilibrated MLS Recovery Technique 
For Error Estimation In Energy Norm 

The discretization error is defined as the difference between the exact solution u and 
the finite element solution u^: e = u — u^. Since the exact solution is in practice 
unknown, in general, the exact error can only be estimated. To obtain an estimation 
of e, measures in terms of the energy norm are normally used. The Zienkiewicz-Zhu 
error estimator defined as 

l|e|| 2 ~ ||e es || 2 = / (v* - <r h ) T B- 1 (v* - <r h ) dQ (12) 

relies on the recovery of an improved stress field <r*, which is supposed to be more 
accurate than the FE solution cr h , to obtain an estimation of the error in energy 
norm ||e e5 ||. 

The proposed recovery technique is based on a moving least squares proce- 
dure which provides a continuous recovered stress field [27] and will be denoted 
as MLSCX. The satisfaction of boundary and internal equilibrium has been con- 
sidered in the formulation, aiming to create a statically admissible stress field that 
would provide accurate error estimates in the FEM and XFEM frameworks. 
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3.1 MLS recovery 

The MLS technique is based on a weighted least squares formulation biased towards 
the test point where the value of the function is asked. The technique considers a 
polynomial expansion for each one of the components of the recovered stress field in 
the form: 

^*( x ) = P(x)a*(x) i = xx,yy,xy (13) 
where p represents a polynomial basis and a are unknown coefficients 

p(x) = {1 x y x 2 xy y 2 . . .} (14) 
a;(x) = {a 0t (x) ai 4 (x) a 2j (x) a 3j (x) a 4t (x) a 5t (x) . . .} T (15) 
For 2D, the expression to evaluate the recovered stress field reads: 



<r*(x) = { a* yy (x) > = P(x)A(x) 
.<( x ). 



'p(x) 
p(x) 
p(x) 



xx (x) 

a yy (x) } (16) 



The format of (16), considering the three components of the stress vector in a 
single equation, will result useful to impose the constraints required to satisfy the 
equilibrium equations. 

Suppose that % is a point within f] x , being Q x the support corresponding to a 
point x defined by a distance (radius) Rq x . The MLS approximation for each stress 
component at % is given by 

<r*(x, X) = p(x)a*(x) V% GO x , i = xx, yy, xy (17) 

To obtain the coefficients A we have adopted the Continuous Moving Least 
Squares Approximation described in [27]. The following functional is minimised: 



J(x) = / W(x- X ) W* (x, x) ~ ° h (X)] 2 dX 



(18) 



Evaluating dJ/dA = results in the linear system M(x)A(x) = G(x) used to 
evaluate A, where 



G 



(x)= f W (x - X ) P T ix) P (X) dX 
(x) = f W (x - X ) P T ix) <r h ix) d X 



(19) 



Assuming that there are n sampling points of coordinates Xi (I — l---^) within 
the support of x, with weight Hi and being |J(xz)| the jacobian determinant, the 
expressions in (18, 19) can be numerically evaluated as 

n 

j{*) = Y J wi*-xiWi^xi)-<T\xi)] 2 \iixim 

1=1 

n 

M (x) = J2 W (x - x,) P T ( Xl ) P ixi) \JiXiW (20) 

1=1 

G (x) = W (x - xi) P T (X,) ^ iXi) \J(Xi)\Hi 



The integration points for the numerical evaluation of the integrals in the above 
equations correspond to the integration points within f] x used in the FE analysis, for 
which the stress field is already available. In (19) W is the MLS weighting function, 
which in this paper has been taken as the fourth-order spline, commonly used in the 
MLS related literature: 

Jl_ 65 2 + 85 3_ 35 4 if | s |<l 

w(x -* ) = \o if| s |>i (21) 

where s denotes the normalised distance function given by 

s = (22) 

In the more commonly used Discrete MLS approach [27] the functional J(x) 
would be defined as: 



J (x) = £ W (x - xi) y (x, xi) - <r h (X,)] 2 (23) 

1 = 1 

This approach would thus produce similar expressions to the equations shown 
in (20) with the difference that, in the continuous approach, each of the sampling 
points Xi i s weighted by its associated area \J(xi)\Hi- O ur numerical experience has 
shown that the Continuous MLS approximation used in this paper is more accurate 
than the discrete approximation, especially when the distribution of the sampling 
points is not uniform within the support. 

Continuity in <r* is directly provided by the MLS procedure previously described 
because the weighting function W ensures that stress sampling points leave or enter 
the support domain in a gradual and smooth manner when x moves [27]. The 
following sections are devoted to the satisfaction of the equilibrium equations. 

3.2 Satisfaction of the boundary equilibrium equation 

The boundary equilibrium equation must be satisfied at each point along the con- 
tour. In [23, 16, 8], where an SPR-based technique was used, the authors enforced 
the satisfaction of the boundary conditions in patches along the boundary using 
Lagrange Multipliers to impose the appropriate constraints between the unknown 
coefficients to be evaluated. However, this approach produces discontinuities in a 
MLS formulation as we move from a support fully in the interior of the domain to 
a support intersecting the boundary. 

In order to avoid the introduction of discontinuities in the recovered field, we 
have followed a nearest point approach that introduces the exact satisfaction of the 
boundary equilibrium equation in a smooth continuous manner. As the constraint 
is smoothly introduced there is no jump when the support does not longer intersects 
T. For a point x e Q whose support Q x intersects the boundary T, the equilibrium 
constraints are considered only in the closest points Xj £ T on the boundaries within 
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the support of x, as shown in Figure 3. Note that we can have more than one nearest 
point for a given support, as is the case for a point x approaching a corner where 
we take one point for each side of the corner (see Figure 3). In this case, two 
different points have to be considered on the boundary to avoid jumps induced by 
the different boundary conditions when crossing the diagonal that bisects the corner. 




Figure 3: MLS support with boundary conditions applied on the nearest boundary 
points. 

Let us express the stress vector cr*(x, %) in a coordinate system xy aligned with 
the contour at Xj suc h that x is the outward normal vector, rotated an angle a with 
respect to x: 

flr*(x,x)=R(a)<T*(x,x) (24) 
where R is the stress rotation matrix 





^XX 




cos 2 a 


sin 2 a 


sin(2a) 




R = 


r yy 




sin 2 a 


cos 2 a 


— sin(2a) 


(25) 




_ r xy_ 




-sin(2a)/2 


sin(2a)/2 


cos(2a) 





The MLS functional expressed in its continuous version and incorporating the 
boundary constraints reads: 

n 

J(x) = W ( x - Xi) [°* (x, Xi) " a* (*,)] 2 \J(Xi)\Hi+ 

1=1 

nbc 

E^( x -^)[^(x,Xi)-af fe)] 2 (26) 
i=i 

= E W ( x - Xi) [P (X) A (x) - a h ( Xl )} 2 1 J( Xl ) \H l+ 

1 = 1 

nbc 

J]^( x - Xj) h(«)P(Xi)A (x) - of (xj)] 2 i = xx, xy 

3=1 
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where nbc is the number of points % • on the boundary where the known bound- 
ary constraints a? x (in general, those would be the normal o%x and tangential 
&xy stresses) are considered. Evaluating dJ/dA = results in the linear system 
M(x)A(x) = G(x) used to evaluate A, where, in this case 



M = J2 W ( x - Xi) P T (Xi) P (Xi) \J(Xi)\Hi+ 

1=1 

nbc 

^2w(^- Xj )P T ( Xj )rTr,P( Xj ) (27) 

3=1 

n 

G = J2 W (x - xi) P T {Xi) ° h (Xi) \J(Xi)\Hi+ 

1 = 1 

nbc 

£ " X 3 )^ T {xM of( Xj ) (28) 
In the previous equations W is a weighting function defined as: 

ur( ^ ~ Xj) /--6^ + 8 5 2 -3 5 3 if |s| < 1 

V^(x-x-) = ^ =< ^ 5 ( 29 ) 

5 [0 if|s|>l 

This function has two main characteristics: 

1. W includes the weighting function W such that the term for the boundary 
constraint is introduced smoothly into the functional J(x). As a result, the 
recovered stress field will be continuous in Q 



2. W also includes s -1 such that the weight of the boundary constraint in J(x) 
increases as we approach the boundary (when x Xj s ~^ 0)? therefore <r* 
will tend to exactly satisfy boundary equilibrium as x — >• Xj ( see Figure 4). 
Note that to estimate the error using the numerical integration in (12), the 
value of <r* is never evaluated on the boundary (where s = 0) because the 
integration points considered are always inside the elements. 



3.3 Satisfaction of the internal equilibrium equation. 

In addition to the enforcement of boundary equilibrium, we will also consider the 
satisfaction of the internal equilibrium equation using the Lagrange Multipliers tech- 
nique. Thus, we will try to enforce the recovered stress field <r* to satisfy the internal 
equilibrium equation 

V • <r* + b = (30) 
The spatial derivatives of <r*, considering (16), are expressed as 

V • <t* = (V • P) A + P (V • A) (31) 
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Figure 4: Satisfaction of boundary equilibrium, a* (x^,%) and a^x^x) are the 
values of cr* (x, %), projected along the direction normal to boundary T at /, in 
the supports Qa and Qb of the points A and 5, whose nearest point on V is /. t n 
represents the normal tractions applied on V. Note that a* (x, Xi) tn ( x /) although 
(j* (x B ,x 7 ) is more accurate than a*(x^,x 7 ). Thus, as x 4 X/, a* (x, % 7 ) — >> 
t n (xj) and, similarly the value of the stresses evaluated at the center of the support 
a* (x, x) t n (xj) 



The first terms in (31) can be directly evaluated differentiating the polynomial 
basis. Previous works [19, 20, 21] have only considered the first term in the sat- 
isfaction of the appropriate equations, thus only providing a pseudo-satisfaction of 
these equations [21]. Therefore, the second term in (31) must also be obtained. To 
evaluate it, we differentiate the linear system MA = G: 

(V • M) A + M (V • A) = V • G (32) 

Evaluating V • A from (32), replacing in (31) and expanding leads to: 

^.(* PM-f)A + PM-f.E lA+f , ,33) 
OX \ox ox J ox 

d °" = - PM- A + PM-f = E,„A + f, , (34) 



dy \dy dy J dy 

where the partial derivatives of M and G with respect, for example, to x are 

i=i 

^d_m^A pT(Xj)rf[;PiXj) (35) 

.7=1 

1=1 

^9W(^pA pT(Xj)r r ar(Xj) (36) 
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where, differentiating (21, 29), 

dW (x - x) dW (x - x) ds 



dx ds dx 



9^ (x - Xj ) = dW (x - Xj ) ds 
dx ds dx 

In these equations ds/dx can be obtained from (22) or, alternatively, from (43) 
for the case shown in the next section. Equations (33, 34) are expressed as a function 
of A, so, we can write the two terms of the internal equilibrium equation (30) as a 
function of the vector of unknowns A: 

(J^xx,x H~ ~^xy,y) A ~l~ (fxx,x H~ fxy,y) ^x (39) 

(J^xx,y + ^yy,y) A + {fxy,x + fyy,y) ~\~ b y = (40) 

where D^- (i = xx,yy,xy and j = x,y) represents the row in \3j corresponding 
to the i th component of the stresses. These expressions define the constraints be- 
tween the coefficients A required to satisfy the internal equilibrium equation at x. 
Lagrange Multipliers are used to impose these constraint equations. 

The use of the Lagrange Multipliers technique to impose the equilibrium con- 
straint (39, 40) in (26) leads to the following system of equations: 



"M 


C T " 




A 




G 


C 







A 




D 



(41) 



where C and D are the terms used to impose the constraint equations and A is the 
vector of Lagrange Multipliers. 

However, in (32) it was assumed that A is evaluated solving MA = G, although, 
operating by blocks in (41) the following system of equations is obtained: 

MA + C T A = G (42) 

Hence, in the formulation proposed in this paper we have neglected the term C T A 
when evaluating the partial derivatives of A. Evidently, this implies that the internal 
equilibrium equation is not fully satisfied, leading to a nearly exact satisfaction of 
the internal equilibrium equation. As described in the numerical examples, this 
approximation represents an enhancement with respect to the pseudo satisfaction 
of equilibrium [21]. 

References [23, 28] show that the error estimator in (12) would produce an upper 
error bound if <r* is statically admissible. The MLSCX recovery technique produces 
a continuous stress field where the internal equilibrium equation is not fully satisfied. 
Hence <r* is continuous and nearly equilibrated and, thus, nearly statically admis- 
sible. Therefore, although the error estimate provided by the proposed recovery 
technique is very sharp, it is not a guaranteed upper error bound. 



da: 



dx 
da, 



m +^k + b x = 



xy 



dx 



+ 



dy 

da: 



yy 



dy 
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3.4 Visibility 



For problems with re-entrant corners a visibility criterion is used to modify the 
normalised distance s in (22). The standard weight function depends on the distance 
between the central point of the support and the sampling points, decreasing as the 
sampling points are located farther from the centre [5]. 

Consider a domain with a re-entrant corner as shown in Figure 5. The value 
of the weight function for a sampling point Xu considering a centre point x whose 
support contains the singularity at % A , diminishes with the visibility of Xi from x 
such that, for points that cannot be directly viewed from x, instead of (22), the 
following equation is used 




Figure 5: Domain with re-entrant corner. 



3.5 Stress splitting for singular problems. 

It is well known that smoothing techniques perform badly when the solution contains 
a singularity. In [7, 29] a technique that decomposes the stress field in singular and 
smooth parts in order to improve the accuracy of SPR-based error estimators was 
proposed. The authors indicated that the exact stress field a corresponding to a 
singular problem can be expressed as the contribution of a smooth stress field, cr srn0) 
and a singular stress field, o- sing 

& & smo & sing (^^) 

Hence, the recovered stress field for this kind of problems can be expressed as 
the contribution of a smooth and a singular recovered stress fields 

°* = <r*8mo + a *8ing ( 45 ) 

To obtain an accurate approximation of the singular part we use the interaction 
integral, as shown in [30], to compute a good estimation of the GSIFs K\ and K\\. 
Then, using the estimated values K{ and K{ Y we can evaluate a singular recovered 
stress field cr* sing from (8). 

Assuming that (T* si is a good approximation of the singular part (T S i ng , a FE- 
type representation of the smooth part cr h srno is given by 

^smo = v h - <r* sing (46) 
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In [7, 29] an SPR-based recovery technique was used to smooth the discontin- 
uous stress field cr h srno . In this paper, we use the moving least squares procedure 
previously described to recover the smooth part of the solution cr* mo . In [7, 29] 
the stress splitting procedure was only used in a small area around the crack tip. 
In the procedure proposed herein the stress splitting is used in the whole domain 
of the problem in order to avoid discontinuities along the blending zone. Thus, the 
boundary tractions to be considered for the satisfaction of the boundary equilibrium 
equation in the smooth problem are: 

t S mo = t — t sing (47) 

where t* si are the projection of (T* si . It must be taken into account that the 
crack faces are treated as any other Neumann boundary where satisfaction of the 
boundary equilibrium equation will be imposed. 

Note that (T* ing is equilibrated and continuous, therefore, the resulting recovered 
stress field <r* = cr* sing + (T* mo only has small lacks of internal equilibrium in <r* mo 
induced by the recovery process. 



3.6 Adaptive strategy 

The refinement of the mesh using the error estimate as the guiding parameter con- 
siders an stopping criterion that checks the value of the estimated error against a 
prescribed or desired error. If the estimated error is higher than the desired error 
then the mesh is refined. Several procedures to perform the refinement are available 
in the literature. To define the size of the elements in the new mesh we follow the 
adaptive process described in [31, 32, 33] which minimises the number of elements 
in the new mesh. This criterion is equivalent to the traditional approach of equally 
distributing the error in each element of the new mesh as proven in [34, 35]. 



4 Numerical Examples 

In this section numerical tests using 2D benchmark problems with exact solution 
are used to investigate the quality of the proposed error estimation technique. The 
first three problems (smooth and singular) consider a FEM approximation whilst 
the fourth problem is solved using an XFEM formulation. For all the models we 
assume a plane strain condition. Sequences of meshes with linear (TRI3), quadratic 
(TRI6) triangles and linear (QUAD4), quadratic (QUAD8) quadrilaterals elements 
are considered for the analyses. Uniform and /i-adaptive refinements have been 
used. The /i-adaptive refinement is based on element splitting using multipoint 
constraints (MPC) to impose C° continuity at hanging nodes. Quadrature rules of 
1, 3, 2 x 2 and 3 x 3 Gauss points are used for TRI3, TRI6, QUAD4 and QUAD8 
elements, respectively. A support size with a radius two times the average size 
of the surrounding elements is used to perform the MLS recovery. 19 sampling 
points in triangular elements and 25 sampling points in quadrilaterals are used for 
an accurate numerical evaluation of (12) in order to avoid the effect of numerical 
errors due to integration. The computational cost of the proposed technique could 
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be alleviated by evaluating (12) using quadrature rules with fewer integration points 
at the expense of introducing errors due to integration in the procedure. The MLS 
basis functions used in the recovery are polynomials p one order higher than the 
corresponding FE displacement basis. 

The performance of the technique is evaluated using the effectivity index of the 
error in energy norm, both at global and local levels. Globally, we consider the value 
of the effectivity index 9 given by 

9 = ^4 (48) 
||e|| ' f 

where ||e|| denotes the exact error in energy norm, and ||e e J| represents the evaluated 
error estimate. At element level, the distribution of the local effectivity index D, its 
mean value m(\D\) and standard deviation cr(D) is analysed, as described in [7]: 

d = e e - i if e e > i \\ e e II 

0-1-1 if «•<! With tr = Wl (49) 

Qe II II 

where superscript e denotes evaluation at element level. 

The h-adaptive refinement procedure considering the error in quantities of inter- 
est is implemented based on previous adaptive procedures using the error in energy 
norm. The technique aims to minimise the number of elements to get the target 
error by equally distributing the element error in the mesh. 

4.1 2x2 square with a 3rd-order polynomial solution 

The 2x2 square model shown in Figure 6 is analysed, with material parameters 
E = 1000 for the Young's modulus and v = 0.3 for the Poisson's ratio. Dirichlet 
boundary conditions are indicated in the figure. The problem is defined such that 
the exact displacement solution is given by 

u(x, y) = x + x 2 — 2xy + x 3 — 3xy 2 + x 2 y (50) 
v(x, y) = -y- 2xy + y 2 - 3x 2 y + y 3 - xy 2 (51) 

The exact values of the stress components are applied along the Neumann bound- 
ary denoted by a dashed line in Figure 6. These stresses can be derived from the 
exact displacement field under plane strain condition, and read 

E 

a xx = Y^t 1 + 2x ~ 2 V + 3 ^ 2 - 3y 2 + 2xy) (52) 
E 

Gyy = ~ 2X + ^ ~ + ~ ( 53 ) 

The following body forces must be applied to satisfy equilibrium: 

E 

b x (x,y) = -^-^(l + y) (55) 
E 
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Figure 6: 2x2 square plate. 

We have used this problem to analyse the influence of different implementations 
of the MLS recovery technique in the error estimate, considering the following cases: 



• MLS: Plain Moving Least Squares recovery 



• MLS+BE: MLS technique with the boundary equilibrium enhancement de- 
scribed in Section 3.2 



• MLS+BE+PIE: MLS technique with the boundary equilibrium enhance- 
ment described in Section 3.2 and the pseudo satisfaction of the internal equi- 
librium equation 



• MLSCX: Technique proposed in this paper 

The results for the plain MLS case will be used as reference. The other three cases 
represent implementations which increasingly approach the full satisfaction of the 
equilibrium equations. Figures 7 to 10 show the effect ivity of the error estimation 
vs. the number of degrees of freedom (dof) using these four implementations for 
/i-adaptive meshes. 

These figures clearly show that the satisfaction of boundary equilibrium (curves 
MLS+BE) plays the most important role towards the enforcement of equilibrium 
and, therefore, an improvement on the accuracy of the error estimator when com- 
pared with the MLS curve. The additional pseudo-satisfaction of internal equi- 
librium (curves MLS+BE+PEI) does not improve, and sometimes provides worse 
effectivities than boundary equilibrium constraints, as it can be seen in Figure 10. 
From Figures 7 to 10 we can see an increase in the accuracy for the MLSCX curves 
with respect to the other curves, with effectivities very close to 6 = 1. 

Figure 11 shows the evolution with respect to mesh refinement of the global 
effectivity index the mean absolute value m(\D\) and standard deviation cr(D) 
of the local effectivity index for the different types of elements considered. Note 
that with the proposed technique we obtain very accurate values of 9 and the error 
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Figure 7: 2 x 2 square with /^-adaptive meshes and TRI3 elements. Evolution of 9 
for different recoveries 



0.99 



0.98 - 



0.97 




-MLS 
- MLS+BE 
-MLS+BE+PIE 
MLSCX 



10 s 



10 4 



MLS 



MLS+BE 



dof 



dof 



366 0.966 366 0.993 

1076 0.977 1068 0.987 

3120 0.988 3068 0.993 

8204 0.992 7920 0.996 



MLS+BE+PIE 



MLSCX 



dof 



dof 



366 0.992 366 0.996 

1068 0.987 1068 0.990 

3038 0.992 3038 0.994 

7866 0.996 7956 0.997 



dof 



Figure 8: 2 x 2 square with /i-adaptive meshes and TRI6 elements. Evolution of 9 
for different recoveries. 



estimate converges to the exact value with the increase of the number of degrees of 
freedom. For this example, the best results are obtained with quadratic elements. 

Figure 12 shows the distribution of the local effectivity on a set of TRI3 meshes, 
Figure 13 displays the same results for a set of QUAD4 meshes. In both cases we can 
observe a quite homogeneous distribution of the local effectivity inside the domain 
and good results along the boundary of the problem. In addition, D decreases for 
finer meshes and we always have values within a very narrow range. Figures 14 and 
15 show a similar behaviour for quadratic elements. 



4.2 Thick- wall cylinder subjected to an internal pressure. 

The geometrical model for this problem is shown in Figure 16. Due to symmetry 
conditions, only one part of the section is modelled. 

The exact solution for this problem is given by the following expressions. For a 
point (x, t/), c — b/a : r = \J x 1 + y 2 the radial displacement is given by 
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Figure 9: 2 x 2 square with /i-adaptive meshes and QUAD4 elements. Evolution of 
for different recoveries 
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Figure 10: 2x2 square with /i-adaptive meshes and QUAD8 elements. Evolution 
of 6 for different recoveries. 
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Figure 11: 2x2 square with uniformly refined meshes. Evolution of the effectivity 
index 9 for different element types. 



Stresses in cylindrical coordinates are 
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Figure 12: 2x2 square with TRI3 elements. Distribution of the effect ivity index D 
in uniformly refined meshes. 
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Figure 13: 2x2 square with QUAD4 elements. Distribution of the effect ivity index 
D in uniformly refined meshes. 
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Figure 14: 2x2 square with TRI6 elements. Distribution of the effectivity index D 
in uniformly refined meshes. 



° r c 2 -ll r 2 / n P 



2 \ a z = 2v- (58) 



1 c 2 — 1 \ r 2 

Figure 17 shows the effectivity values obtained with the MLSCX in the thick- wall 
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Figure 15: 2x2 square with QUAD8 elements. Distribution of the effect ivity index 
D in uniformly refined meshes. 
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Figure 16: Thick- wall cylinder subjected to internal pressure. 

cylinder using uniformly refined meshes. The results obtained are similar to those 
previously shown for the square plate. Compensations between underestimated 
and overestimated areas of the domain might result in misleading values of the 
global value 9. To take this into account we consider the parameters m{\D\) and 
(j{D) which are expected to decrease when we increase the level of refinement. In 
Figure 17 m(\D\) and cr(D) decrease when we increase the refinement showing a 
good performance of the error estimator. The results for this problem show that the 
proposed technique provides an accurate estimate of the exact error in energy norm. 
On the other hand, for this example the best values are obtained when considering 
linear elements, thus, there is not a strong correlation between the order of the 
approximation and the performance of the error estimator. 

Figure 18 shows the distribution of the exact error in energy norm for the same 
meshes. The higher errors are located close to the inner radius of the cylinder. The 
error decreases as we increase the number of degrees of freedom. 

Figure 19 shows the distribution of the local effectivity index Dina sequence of 
meshes with linear triangular elements. The figure displays a quite uniform distribu- 
tion of the local effectivity at each mesh, within the range between [—0.54, 0.43] for 
the coarsest mesh. It is worth noting that the local values improve as we refine the 
meshes and that, for the last mesh in the sequence, the local effectivity is now within 
the range [—0.26,0.17]. Some radial patterns in the distribution of D can be seen 

20 




Figure 17: Thick-wall cylinder with uniformly refined meshes. Global indicators 0, 
m(\D\) and a(D). 
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Figure 18: Thick- wall cylinder with TRI3 elements. Distribution of the exact error 
in energy norm in uniformly refined meshes. 



as we increase the number of dof, which are attributed to local mesh configurations. 

Figure 20 shows the distribution of the local effectivity index D in a sequence 
of meshes with linear quadrilateral elements (QUAD4). Again, we can see a quite 
uniform distribution of the local effectivity at each mesh and that the local effectivity 
improves as the mesh is refined. This behaviour indicates that the proposed error 
estimator performs nicely at element level, which is important when guiding adaptive 
processes. 
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Figure 19: Thick-wall cylinder with TRI3 elements. Distribution of the effectivity 
index D in uniformly refined meshes. 
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Figure 20: Thick- wall cylinder with QUAD4 elements. Distribution of the effectivity 
index D in uniformly refined meshes. 

4.2.1 Influence of support size. 

One of the parameters that affects the performance of the proposed error estimator is 
the radius Rq x that defines the MLS support at a point x, also known as the domain 
of influence. The idea is to define a support with a number of sampling points large 
enough to be able to solve the MLS fitting and to obtain an accurate polynomial 
expansion of the stresses, but not too large that we risk excessively smoothing the 
stress field and no longer describing the local behaviour of the solution. Moreover, 
larger supports means more computational effort as more sampling points should be 
considered. 

In order to fix the domain of influence at a particular point we first evaluate the 
average size of the elements surrounding each node of the mesh. Then, we define 
the radius of the support at nodes as i?^ x (x^) = fc/(x^) where k is a constant that 
takes positive values and /(x^) is the average size of the elements containing node i. 
Once the value of Rq x is evaluated at nodes, the value of Rq x at any point x within 
an element is interpolated from the nodes using the displacement shape functions. 
Note that as Rq x is a function of x. Its definition has been used in the derivatives 
of s defined in (22) (or alternatively in (43)) required for the evaluation of (37, 38). 

Figure 21 shows the global results for a sequence of TRI3 elements considering 
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different values of k. Figure 22 shows the same results for QUAD4 meshes. Note 
that for small supports the values of the local indicators are less accurate even if the 
values for the global indicator are closer to one. A good balance between accuracy 
and local definition of the smoothing function is obtained for k = 2, which is the 
value considered in the examples presented herein. 




dof dof 

Figure 21: Thick- wall cylinder with uniformly refined TRI3 meshes. Global indica- 
tors m(\D\) and &(D) for different values of k. 



4.3 Westergaard problem — FEM solution. 

To evaluate the performance of the proposed technique for singular problems we 
consider the Westergaard problem [7, 36] as it has an exact analytical solution. The 
Westergaard problem corresponds to an infinite plate loaded at infinity with biaxial 
tractions a xoo = a yoo = and shear traction r^, presenting a crack of length 2a as 
shown in Figure 23. Combining the externally applied loads we can obtain different 
loading conditions: pure mode I, pure mode II or mixed mode. 

The numerical model corresponds to a finite portion of the domain (a = 1 and 
b = 4 in Figure 23). The applied projected stresses for mode I are evaluated from 
the analytical Westergaard solution [36] : 
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Figure 22: Thick- wall cylinder with uniformly refined QUAD4 meshes. Global indi- 
cators 0, m(\D\) and cr(D) for different values of k. 
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and for mode II: 
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where the stress fields are expressed as a function of x and ?/, with origin at the 
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Figure 23: Westergaard problem. Infinite plate with a crack of length 2a under 
uniform tractions (biaxial) and r^. Finite portion of the domain fi , modelled 
with FE. 

centre of the crack. The parameters i, m, n and (j) are defined as 

t — {x + iy) 2 — a 2 = (x 2 — y 2 — a 2 ) + i(2xy) = m + in 
m = Re(t) = Re(z 2 — a 2 ) = x 2 — y 2 — a 2 
n = Im(t) = (z 2 — a 2 ) = 2xy 

Arg(f) Arg(m — m) with (j) G [— 7r, 7r] , i 2 = — 1 
For the problem analysed, the exact value of the SIF is given by 

Kl,ex &ocV™ Kll,ex T'ocV™ (62) 

Material parameters are Young's modulus E = 10 7 and Poisson's ratio v = 0.333. 
We consider loading conditions in pure mode I with = 100 and = 0, pure 
mode II with and 100, and mixed mode with 100 and 100. 

To evaluate the SIF needed for recovering the singular part we use an equivalent 
domain integral technique, with a plateau function with radius r q = 0.9 for the 
extraction [37]. 

Figure 24 shows the evolution with respect to mesh refinement of the global 
parameters 6 ) m(\D\) and cr(D) for different element types. In the figure, the global 
effectivity converges to the theoretical value of 9 = 1 and both m(\D\) and cr(D) 
decrease with an increase of the number of dof. The performance of the proposed 
technique indicates an accurate error estimation for the meshes analysed. 

Figures 25 and 26 show the distribution of the local effectivity index D in a 
sequence of TRI3 and QUAD4 meshes respectively. The splitting of the stress field 
into singular and smooth parts helps to recover a highly accurate stress field in 
the vicinity of the singular point. The distribution of the local effectivity index is 
homogeneous within the mesh and the values for D decrease as we refine. 

Because local error estimation techniques cannot take into account the pollution 
error due to the singularity, we can notice areas of the domain where the error is 
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Figure 24: Westergaard problem under mode I with FEM /i-adapted meshes. Global 
indicators 0, m(\D\) and cr(D). 



underestimated in this example, especially in the first meshes. The effect of pollution 
error is partially overcome by the use of //.-adaptive refinement (or enriched meshes 
as it is shown in the next section). 



4.4 Westergaard problem - XFEM solution. 

Let us now consider the Westergaard problem from the previous section, solved 
using an enriched finite element approximation. In the numerical analyses, we use 
a geometrical enrichment defined by a circular fixed enrichment area £>(x , r e ) with 
radius r e = 0.5, with its centre at the crack tip x as proposed in [38]. For the 
extraction of the SIF we define a plateau function with radius r q = 0.9 as in the 
FEM case. Bilinear elements are considered in the models. For the numerical 
integration of standard elements we use a 2 x 2 Gaussian quadrature rule. The 
elements intersected by the crack are split into triangular integration subdomains 
that do not contain the crack. Alternatives which do not require this subdivision 
are proposed in [39, 40]. We use 7 Gauss points in each triangular subdomain, and 
a 5 x 5 quasipolar integration in the subdomains of the element containing the crack 
tip [38], see Figure 2. We do not consider correction for blending elements. Methods 
to address blending errors are proposed in [41, 42, 43, 44]. 

Figure 27 shows the evolution with respect to mesh refinement of the global 
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Figure 25: Westergaard problem under mode I with FEM and /^-adapted meshes of 
TRI3. Distribution of the effectivity index D. 
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Figure 26: Westergaard problem under mode I with FEM and /^-adapted meshes of 
QUAD4. Distribution of the effectivity index D. 



parameters 0, m{\D\) and a(D) for the structured meshes of enriched QUAD4 ele- 
ments. The curves represent the values obtained for the Westergaard problem under 
mode I, mode II and mixed mode loading conditions. In the figure, the global ef- 
fectivity converges to the theoretical value of 9 = 1 and both m(\D\) and cr(D) 
decrease with an increase of the number of dof. The results show that the proposed 
technique provides a sharp estimate of the true error. 

Figure 28 shows the distribution of D in the second mesh (1895 dof) of the 
sequence of structured meshes for all the three loading modes. The results indicate 
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Figure 27: Westergaard problem with XFEM and structured meshes of QUAD4. 
Global indicators 0, m(\D\) and cr(D). 

a quite uniform distribution of the local effectivity. The values of D indicate that 
the error at element level is accurately evaluated even where standard recovery 
techniques would produce the worst results: along the Neumann boundary, the 
crack faces and around the crack tip. 

For the case of non structured meshes the results for the same global parameters 
previously considered are shown in Figure 29. The local effectivity at element level 
for this meshes is depicted in Figure 30. There is a similar behaviour to that seen 
for structured meshes. In general, the proposed technique exhibits an excellent 
performance when used to estimate the error in the XFEM approximations analysed. 

In Figure 31 we compare the results of the MLSCX with those of the SPRCX 
recovery procedure. In this case both techniques give values in the same order of 
magnitude. 



5 Conclusions 

In this paper, the use of an equilibrated moving least squares recovery technique for 
FEM and XFEM problems has been investigated. The proposed technique uses a 
MLS approach to provide a continuous recovered stress field that enforces boundary 
equilibrium constraints. It also imposes a very accurate satisfaction, although not 
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Figure 28: Westergaard problem with XFEM and structured meshes of QUAD4: a) 
mode I, b) mode II and c) mixed mode. Distribution of the effectivity index D. 

fully exact, of the internal equilibrium equation. Moreover, for singular problems it 
decomposes the stress field into two different parts, singular and smooth, in order to 
enable the technique to describe the singular behaviour of the solution. A visibility 
criterion is used near reentrant corners and cracks to properly define the weight of 
the sampling points within the support. 

The technique presented here has been validated using four different examples 
with known analytical solution. The numerical results have shown the accuracy of 
the proposed technique, which provides values of the effectivity index that converge 
and are very close to the theoretical value 6—1. The distribution of the local 
effectivity at the elements is homogeneous for the tests considered, and the mean 
value m(\D\) and standard deviation a(D) decrease as we increase the number of 
dof. The obtained MLS recovered field is not fully statically admissible, thus, the 
procedure does not guarantee the upper bound property. For this reason, it nearly 
bounds the exact error but not always yields an effectivity index greater than one, 
as clearly seen in the first example. In any case, the numerical results show that for 
the examples presented the proposed technique yields sharp error estimates, which 
are very accurate when compared with previous MLS approaches. Extension of this 
work to 3D problems is feasible given that the SIF along the crack front is evaluated 
with sufficient accuracy. It is known that in 3D problems the evaluation of the SIF 
is less accurate. In [7] the influence of the accuracy in the evaluation of the SIF in 
the error estimator is investigated. 
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Figure 29: Westergaard problem with XFEM and non structured meshes of QUAD4. 
Global indicators m(\D\) and cr(D). 
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Figure 30: Westergaard problem with XFEM and non structured meshes of QUAD4: 
a) mode I, b) mode II and c) mixed mode. Distribution of the effectivity index D. 

knowledged. 

Support from the EPSRC grant EP/G042705/1 "Increased Reliability for Industri- 



30 



1.2 



1.15 



1.1 



1.05 



-e- MLSCX 
-a- SPRCX 



1.2 




10 3 



10 4 
dof 



- 1.15 



1.1 



1.05 



10 5 



- MLSCX 

- SPRCX 




10 s 



10 4 
dof 



10 5 



1.2 



1.15 



1.1 



1.05 



0.95 



-0- MLSCX 
-B- SPRCX 




10 6 



10 4 
dof 



10° 



Figure 31: Westergaard problem with XFEM and structured meshes of QUAD4. 
Effectivity index for the MLSCX and SPRCX recovery techniques. 
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